Załącznik 1
Dwa kody MATLABa obliczające izogeometryczną L2 projekcję bitmapy
(zob. rozdział Implementacja w MATLABie problemu projekcji bitmapy )
KOD 1
function bitmap_param(filename,nxx,pxx,nyy,pyy)%, knot_vectorx, knot_vectory)
%funkcja liczaca ilosc funkcji bazowych
compute_nr_basis_functions = @(knot_vector,p) size(knot_vector, 2) - p - 1;
knot_vectorx = simple_knot(nxx,pxx);
knot_vectory = simple_knot(nyy,pyy);
%filename = './2019-06-20 17.49.00.tif';
%to juz jest RGB
X = imread(filename);
R = X(:,:,1);
G = X(:,:,2);
B = X(:,:,3);
ix = size(X,1);
iy = size(X,2);
px = compute_p(knot_vectorx);
py = compute_p(knot_vectory);
elementsx = number_of_elements(knot_vectorx,px);
elementsy = number_of_elements(knot_vectory,py);
nx = number_of_dofs(knot_vectorx,px);
ny = number_of_dofs(knot_vectory,py);
A = sparse(nx*ny,nx*ny);
FR = zeros(nx*ny,1);
FG = zeros(nx*ny,1);
FB = zeros(nx*ny,1);
%calki z B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
%(i,k=1,...,Nx; j,l=1,...,Ny)
%petla po elementach w osi x
for ex = 1:elementsx;
% zakres funkcji niezerowych nad elementem
[xl,xh] = dofs_on_element(knot_vectorx,px,ex);
% zakres elementu (lewy i prawy brzeg po x)
[ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
%petla po elementach w osi y
for ey = 1:elementsy
% zakres funkcji niezerowych nad elementem
[yl,yh] = dofs_on_element(knot_vectory,py,ey);
% zakres elementu (lewy i prawy brzeg po x)
[ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
% jakobian - rozmiar elementu
Jx = ex_bound_h - ex_bound_l;
Jy = ey_bound_h - ey_bound_l;
J = Jx*Jy;
% petla po funkcjach niezerowych nad danym elementem
for bi = xl:xh
for bj = yl:yh
for bk = xl:xh
for bl = yl:yh
% punkty kwadratury w osi x nad elementem
qpx = quad_points(ex_bound_l,ex_bound_h,2*px+2*py+1);
% punkty kwadratury w osi y nad elementem
qpy = quad_points(ey_bound_l,ey_bound_h,2*px+2*py+1);
% wagi kwadratury w osi x nad elementem
qwx = quad_weights(ex_bound_l,ex_bound_h,2*px+2*py+1);
% wagi kwadratury w osi y nad elementem
qwy = quad_weights(ey_bound_l,ey_bound_h,2*px+2*py+1);
% petla po punktach kwadratury
for iqx = 1:size(qpx,2)
for iqy = 1:size(qpy,2)
% definicja funkcji ksztaltu
% B^x_k(x)
funk= compute_spline(knot_vectorx,px,bi,qpx(iqx));
% B^y_l(y)
funl= compute_spline(knot_vectory,py,bj,qpy(iqy));
% B^x_i(x)
funi= compute_spline(knot_vectorx,px,bk,qpx(iqx));
% B^y_j(y)
funj= compute_spline(knot_vectory,py,bl,qpy(iqy));
% B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
fun = funi*funj*funk*funl;
% Calki z B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
% (i,k=1,...,Nx; j,l=1,...,Ny)
int = fun*qwx(iqx)*qwy(iqy)*J;
if (int!=0)
A((bj-1)*nx+bi,(bl-1)*nx+bk) = A((bj-1)*nx+bi,(bl-1)*nx+bk) + int;
endif
endfor
endfor
endfor
endfor
endfor
endfor
endfor
endfor
%Calki BITMAPA(x,y) B^x_k(x) B^y_l(y)
%petla po elementach w osi x
for ex = 1:elementsx;
% zakres funkcji niezerowych nad elementem
[xl,xh] = dofs_on_element(knot_vectorx,px,ex);
% zakres elementu (lewy i prawy brzeg po x)
[ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
%petla po elementach w osi y
for ey = 1:elementsy
% zakres funkcji niezerowych nad elementem
[yl,yh] = dofs_on_element(knot_vectory,py,ey);
% zakres elementu (lewy i prawy brzeg po x)
[ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
%jakobian - rozmiar elementu
Jx = ex_bound_h - ex_bound_l;
Jy = ey_bound_h - ey_bound_l;
J = Jx * Jy;
%petla po funkcjach niezerowych nad danym elementem
for bi = xl:xh
for bj = yl:yh
%definicja funkcji ksztaltu
% punkty kwadratury w osi x nad elementem
qpx = quad_points(ex_bound_l,ex_bound_h,px*py+1);
% punkty kwadratury w osi y nad elementem
qpy = quad_points(ey_bound_l,ey_bound_h,px*py+1);
% wagi kwadratury w osi x nad elementem
qwx = quad_weights(ex_bound_l,ex_bound_h,px*py+1);
% wagi kwadratury w osi y nad elementem
qwy = quad_weights(ey_bound_l,ey_bound_h,px*py+1);
%petla po punktach kwadratury
for iqx = 1:size(qpx,2)
for iqy = 1:size(qpy,2)
% B^x_k(x)
funk= compute_spline(knot_vectorx,px,bi,qpx(iqx));
% B^y_l(y)
funl = compute_spline(knot_vectory,py,bj,qpy(iqy));
% Calki BITMAPA(x,y) B^x_k(x) B^y_l(y) po kolorach RGB
intR = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(R,qpx(iqx),qpy(iqy));
intG = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(G,qpx(iqx),qpy(iqy));
intB = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(B,qpx(iqx),qpy(iqy));
FR((bj-1)*nx+bi) = FR((bj-1)*nx+bi) + intR;
FG((bj-1)*nx+bi) = FG((bj-1)*nx+bi) + intG;
FB((bj-1)*nx+bi) = FB((bj-1)*nx+bi) + intB;
endfor
endfor
endfor
endfor
endfor
endfor
% rozwiazanie 3 ukladow rownan
RR=A\FR;
GG=A\FG;
BB=A\FB;
%odtworzenie wyniku
%wyzerowanie macierzy odtworzonego wyniku
R1 = zeros(ix,iy);
G1 = zeros(ix,iy);
B1 = zeros(ix,iy);
% petla po funkcjach
for bi = 1:nx
for bj = 1:ny
% petla po pixelach obrazka %na których funkcje sa niezerowe
for i=xx(knot_vectorx(bi)):xx(knot_vectorx(bi+px+1))
% przeliczenie wspolrzednych [1-szerokosc] -> [0-1]
ii = (i-1)/(ix-1);
% B^x_i(x)
funi= compute_spline(knot_vectorx,px,bi,ii);
for j=yy(knot_vectory(bj)):yy(knot_vectory(bj+py+1))
% przeliczenie wspolrzednych [1-wysokosc] -> [0-1]
jj = (j-1)/(iy-1);
% B^y_j(y)
funj= compute_spline(knot_vectory,py,bj,jj);
ff = funi*funj;
R1(i,j) = R1(i,j) + floor(ff*RR((bj-1)*nx+bi));
G1(i,j) = G1(i,j) + floor(ff*GG((bj-1)*nx+bi));
B1(i,j) = B1(i,j) + floor(ff*BB((bj-1)*nx+bi));
endfor
endfor
endfor
endfor
RGB=X;
RGB(:,:,1) = R1;
RGB(:,:,2) = G1;
RGB(:,:,3) = B1;
imshow(RGB);
function resx=xx(x)
resx = floor((ix-1)*x+1);
endfunction
function resy=yy(y)
resy = floor((iy-1)*y+1);
endfunction
function val=bitmp(M,x,y)
val = zeros(size(x));
for i=1:size(x,1)
for j=1:size(x,1)
val(i,j)=M(xx(x(1,i)),yy(y(1,j)));
endfor
endfor
endfunction
%funkcja wyliczajaca stopien wielomianow
function p=compute_p(knot_vector)
%pierwszy wpis w knot_vector
initial = knot_vector(1);
%dlugosc knot_vector
kvsize = size(knot_vector,2);
p = 0;
%sprawdzenie ilosci powtorzen pierwszego wpisu w knot_vector
while (p+2 <= kvsize) && (initial == knot_vector(p+2))
p = p+1;
endwhile
return
endfunction
%funkcja sprawdzajaca poprawnosc knot_vector
function t=check_sanity(knot_vector,p)
initial = knot_vector(1);
kvsize = size(knot_vector,2);
t = true;
counter = 1;
%jesli ilosc powtorzen na poczatku nie jest zgodna ze stopniem wielomianow
for i=1:p+1
if (initial != knot_vector(i))
%zwroc falsz
t = false;
return
endif
endfor
%jesli w srodku knot_vector jest za duzo powtorzen
for i=p+2:kvsize-p-1
if (initial == knot_vector(i))
counter = counter + 1;
if (counter > p)
%zwroc falsz
t = false;
return
endif
else
initial = knot_vector(i);
counter = 1;
endif
endfor
initial = knot_vector(kvsize);
%jesli ilosc powtorzen na poczatku nie jest zgodna ze stopniem wielomianow
for i=kvsize-p:kvsize
if (initial != knot_vector(i))
%zwroc falsz
t = false;
return
endif
endfor
%jesli kolejny element knot_vector jest mniejszy niz poprzedni
for i=1:kvsize-1
if (knot_vector(i)>knot_vector(i+1))
%zwroc falsz
t = false;
return
endif
endfor
return
endfunction
%funkcja wyliczajaca rekurencyjnie funkcje bazowa zgodnie z wzorem Cox-de-Boor
function y=compute_spline(knot_vector,p,nr,x)
%funkcja (x-x_i)/(x_{i-p}-x_i)
fC= @(x,a,b) (x-a)/(b-a);
%funkcja (x_{i+p+1}-x)/(x_{i+p+1}-x_{i+1})
fD= @(x,c,d) (d-x)/(d-c);
%x_i
a = knot_vector(nr);
%x_{i-p}
b = knot_vector(nr+p);
%x_{i+1}
c = knot_vector(nr+1);
%x_{i+p+1}
d = knot_vector(nr+p+1);
%funkcja liniowa dla p=0
if (p==0)
y = 0 .* (x < a) + 1 .* (a <= x & x <= d) .+ 0 .* (x > d);
return
endif
%B_{i,p-1}
lp = compute_spline(knot_vector,p-1,nr,x);
%B_{i+1,p-1}
rp = compute_spline(knot_vector,p-1,nr+1,x);
%(x-x_i)/(x_{i-p)-x_i)*B_{i,p-1}
if (a==b)
%jesli wezly w knot_vector sie powtarzaja trzeba to uwzglednic
y1 = 0 .* (x < a) + 1 .* (a <= x & x <= b) + 0 .* (x > b);
else
y1 = 0 .* (x < a) + fC(x,a,b) .* (a <= x & x <= b) + 0 .* (x > b);
endif
%(x_{i+p+1}-x)/(x_{i+p+1)-x_{i+1})*B_{i+1,p-1}
if (c==d)
%jesli wezly w knot_vector sie powtarzaja trzeba to uwzglednic
y2 = 0 .* (x < c) + 1 .* (c < x & x <= d) + 0 .* (d < x);
else
y2 = 0 .* (x < c) + fD(x,c,d) .* (c < x & x <= d) + 0 .* (d < x);
endif
y = lp .* y1 .+ rp .* y2;
return
endfunction
% Ilosc elementow w knot vector
function n = number_of_elements(knot_vector,p)
initial = knot_vector(1);
kvsize = size(knot_vector,2);
n = 0;
for i=1:kvsize-1
if (knot_vector(i) != initial)
initial = knot_vector(i);
n = n+1;
endif
endfor
endfunction
% Tworzy prosty knot vector bez powotrzen w srodku
function knot = simple_knot(elems, p)
pad = ones(1, p);
knot = [0 * pad, 0:elems, elems * pad];
knot = knot/elems;
endfunction
% Ilosc funkcji bazowych w knot wektorze
function n = number_of_dofs(knot,p)
n = length(knot) - p - 1;
endfunction
function first = first_dof_on_element(knot_vector,p,elem_number)
[l,h] = element_boundary(knot_vector,p,elem_number);
first = lookup(knot_vector, l) - p;
endfunction
function [low,high] = element_boundary(knot_vector,p,elem_number)
initial = knot_vector(1);
kvsize = size(knot_vector,2);
k = 0;
low=0;
high=0;
for i=1:kvsize
if (knot_vector(i) != initial)
initial = knot_vector(i);
k = k+1;
endif
if (k == elem_number)
low = knot_vector(i-1);
high = knot_vector(i);
return;
endif
endfor
endfunction
% Zwraca zakres (indeksy) funckji bedacych niezerowymi na zadanym wektorze wezlow
function [low,high] = dofs_on_element(knot_vector,p,elem_number)
low = first_dof_on_element(knot_vector,p,elem_number);
%poniewaz mamy dokladnie p+1 niezerowych funkcji nad elementem
high = low + p;
endfunction
% Row vector of points of the k-point Gaussian quadrature on [a, b]
function xs = quad_points(a, b, k)
% mapowanie punktow
map = @(x) 0.5 * (a * (1 - x) + b * (x + 1));
switch (k)
case 1
xs = [0];
case 2
xs = [-1/sqrt(3), ...
1/sqrt(3)];
case 3
xs = [-sqrt(3/5), ...
0, ...
sqrt(3/5)];
case 4
xs = [-sqrt((3+2*sqrt(6/5))/7), ...
sqrt((3-2*sqrt(6/5))/7), ...
sqrt((3-2*sqrt(6/5))/7), ...
sqrt((3+2*sqrt(6/5))/7)];
case 5
xs = [-1/3*sqrt(5+2*sqrt(10/7)), ...
-1/3*sqrt(5-2*sqrt(10/7)), ...
0, ...
1/3*sqrt(5-2*sqrt(10/7)), ...
1/3*sqrt(5+2*sqrt(10/7))];
otherwise
xs = [-1/3*sqrt(5+2*sqrt(10/7)), ...
-1/3*sqrt(5-2*sqrt(10/7)), ...
0, ...
1/3*sqrt(5-2*sqrt(10/7)), ...
1/3*sqrt(5+2*sqrt(10/7))];
endswitch
xs = map(xs);
endfunction
% Row vector of weights of the k-point Gaussian quadrature on [a, b]
function ws = quad_weights(a, b, k)
switch (k)
case 1
ws = [2];
case 2
ws = [1, 1];
case 3
ws = [5/9, ...
8/9, ...
5/9];
case 4
ws = [(18-sqrt(30))/36, ...
(18+sqrt(30))/36, ...
(18+sqrt(30))/36, ...
(18-sqrt(30))/36];
case 5
ws = [(322-13.0*sqrt(70))/900, ...
(322+13.0*sqrt(70))/900, ...
128/225, ...
(322+13.0*sqrt(70))/900, ...
(322-13.0*sqrt(70))/900];
otherwise
ws = [(322-13.0*sqrt(70))/900, ...
(322+13.0*sqrt(70))/900, ...
128/225, ...
(322+13.0*sqrt(70))/900, ...
(322-13.0*sqrt(70))/900];
endswitch
% Gaussian quadrature is defined on [-1, 1], we use [a, b]
%ws = ws * (b-a)/2;
endfunction
endfunction
KOD 2
function bitmap_fast(filename,nxx,pxx,nyy,pyy)%, knot_vectorx, knot_vectory)
tic;
%funkcja liczaca ilosc funkcji bazowych
compute_nr_basis_functions = @(knot_vector,p) size(knot_vector, 2) - p - 1;
knot_vectorx = simple_knot(nxx,pxx);
knot_vectory = simple_knot(nyy,pyy);
%to juz jest RGB
X = imread(filename);
R = X(:,:,1);
G = X(:,:,2);
B = X(:,:,3);
ix = size(X,1);
iy = size(X,2);
px = compute_p(knot_vectorx);
py = compute_p(knot_vectory);
elementsx = number_of_elements(knot_vectorx,px);
elementsy = number_of_elements(knot_vectory,py);
nx = number_of_dofs(knot_vectorx,px);
ny = number_of_dofs(knot_vectory,py);
A = sparse(nx*ny,nx*ny);
FR = zeros(nx*ny,1);
FG = zeros(nx*ny,1);
FB = zeros(nx*ny,1);
init=toc
tic;
splinex = zeros(elementsx,nx,px+1);
spliney = zeros(elementsy,ny,py+1);
for ex = 1:elementsx;
% zakres funkcji niezerowych nad elementem
[xl,xh] = dofs_on_element(knot_vectorx,px,ex);
% zakres elementu (lewy i prawy brzeg po x)
[ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
% punkty kwadratury w osi x nad elementem
qpx = quad_points(ex_bound_l,ex_bound_h,px+1);
% wagi kwadratury w osi x nad elementem
qwx = quad_weights(ex_bound_l,ex_bound_h,px+1);
% petla po funkcjach niezerowych nad danym elementem
for bi = xl:xh
% petla po punktach kwadratury
for iqx = 1:size(qpx,2)
splinex(ex,bi,iqx)= compute_spline(knot_vectorx,px,bi,qpx(iqx));
endfor
endfor
endfor
for ey = 1:elementsy;
% zakres funkcji niezerowych nad elementem
[yl,yh] = dofs_on_element(knot_vectory,py,ey);
% zakres elementu (lewy i prawy brzeg po y)
[ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
% punkty kwadratury w osi y nad elementem
qpy = quad_points(ey_bound_l,ey_bound_h,py+1);
% wagi kwadratury w osi y nad elementem
qwy = quad_weights(ey_bound_l,ey_bound_h,py+1);
% petla po funkcjach niezerowych nad danym elementem
for bi = yl:yh
% petla po punktach kwadratury
for iqy = 1:size(qpy,2)
spliney(ey,bi,iqy)= compute_spline(knot_vectory,py,bi,qpy(iqy));
endfor
endfor
endfor
init_splines=toc
tic;
%calki z B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
%(i,k=1,...,Nx; j,l=1,...,Ny)
%petla po elementach w osi x
for ex = 1:elementsx;
% zakres funkcji niezerowych nad elementem
[xl,xh] = dofs_on_element(knot_vectorx,px,ex);
% zakres elementu (lewy i prawy brzeg po x)
[ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
%petla po elementach w osi y
for ey = 1:elementsy
% zakres funkcji niezerowych nad elementem
[yl,yh] = dofs_on_element(knot_vectory,py,ey);
% zakres elementu (lewy i prawy brzeg po x)
[ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
% jakobian - rozmiar elementu
Jx = ex_bound_h - ex_bound_l;
Jy = ey_bound_h - ey_bound_l;
J = Jx*Jy;
% punkty kwadratury w osi x nad elementem
qpx = quad_points(ex_bound_l,ex_bound_h,px+1);
% wagi kwadratury w osi x nad elementem
qwx = quad_weights(ex_bound_l,ex_bound_h,px+1);
% punkty kwadratury w osi y nad elementem
qpy = quad_points(ey_bound_l,ey_bound_h,py+1);
% wagi kwadratury w osi y nad elementem
qwy = quad_weights(ey_bound_l,ey_bound_h,py+1);
% petla po funkcjach niezerowych nad danym elementem
for bi = xl:xh
for bj = yl:yh
for bk = xl:xh
for bl = yl:yh
% petla po punktach kwadratury
for iqx = 1:size(qpx,2)
for iqy = 1:size(qpy,2)
% B^x_k(x)
funk = splinex(ex,bk,iqx);
% B^y_l(y)
funl = spliney(ey,bl,iqy);
% B^x_i(x)
funi = splinex(ex,bi,iqx);
% B^y_j(y)
funj = spliney(ey,bj,iqy);
% B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
fun = funi*funj*funk*funl;
% Calki z B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
% (i,k=1,...,Nx; j,l=1,...,Ny)
int = fun*qwx(iqx)*qwy(iqy)*J;
if (int!=0)
A((bj-1)*nx+bi,(bl-1)*nx+bk) = A((bj-1)*nx+bi,(bl-1)*nx+bk) + int;
endif
endfor
endfor
endfor
endfor
endfor
endfor
endfor
endfor
lhs=toc
tic;
%Calki BITMAPA(x,y) B^x_k(x) B^y_l(y)
%petla po elementach w osi x
for ex = 1:elementsx;
% zakres funkcji niezerowych nad elementem
[xl,xh] = dofs_on_element(knot_vectorx,px,ex);
% zakres elementu (lewy i prawy brzeg po x)
[ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
%petla po elementach w osi y
for ey = 1:elementsy
% zakres funkcji niezerowych nad elementem
[yl,yh] = dofs_on_element(knot_vectory,py,ey);
% zakres elementu (lewy i prawy brzeg po x)
[ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
%jakobian - rozmiar elementu
Jx = ex_bound_h - ex_bound_l;
Jy = ey_bound_h - ey_bound_l;
J = Jx * Jy;
% punkty kwadratury w osi x nad elementem
qpx = quad_points(ex_bound_l,ex_bound_h,px+1);
% punkty kwadratury w osi y nad elementem
qpy = quad_points(ey_bound_l,ey_bound_h,px+1);
% wagi kwadratury w osi x nad elementem
qwx = quad_weights(ex_bound_l,ex_bound_h,py+1);
% wagi kwadratury w osi y nad elementem
qwy = quad_weights(ey_bound_l,ey_bound_h,py+1);
%petla po funkcjach niezerowych nad danym elementem
for bk = xl:xh
for bl = yl:yh
%petla po punktach kwadratury
for iqx = 1:size(qpx,2)
for iqy = 1:size(qpy,2)
% B^x_k(x)
funk = splinex(ex,bk,iqx);
% B^y_l(y)
funl = spliney(ey,bl,iqy);
% Calki BITMAPA(x,y) B^x_k(x) B^y_l(y) po kolorach RGB
intR = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(R,qpx(iqx),qpy(iqy));
intG = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(G,qpx(iqx),qpy(iqy));
intB = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(B,qpx(iqx),qpy(iqy));
FR((bl-1)*nx+bk) = FR((bl-1)*nx+bk) + intR;
FG((bl-1)*nx+bk) = FG((bl-1)*nx+bk) + intG;
FB((bl-1)*nx+bk) = FB((bl-1)*nx+bk) + intB;
endfor
endfor
endfor
endfor
endfor
endfor
rhs=toc
tic;
% rozwiazanie 3 ukladow rownan
RR=A\FR;
GG=A\FG;
BB=A\FB;
%odtworzenie wyniku
%wyzerowanie macierzy odtworzonego wyniku
R1 = zeros(ix,iy);
G1 = zeros(ix,iy);
B1 = zeros(ix,iy);
factor=toc
tic;
funx_tab = zeros(nx,ix);
funy_tab = zeros(ny,iy);
% petla po funkcjach
for bi = 1:nx
% petla po pixelach obrazka %na których funkcje sa niezerowe
for i=xx(knot_vectorx(bi)):xx(knot_vectorx(bi+px+1))
% przeliczenie wspolrzednych [1-szerokosc] -> [0-1]
ii = (i-1)/(ix-1);
% B^x_i(x)
funx_tab(bi,i) = compute_spline(knot_vectorx,px,bi,ii);
endfor
endfor
% petla po funkcjach
for bj = 1:ny
% petla po pixelach obrazka %na których funkcje sa niezerowe
for j=yy(knot_vectory(bj)):yy(knot_vectory(bj+py+1))
% przeliczenie wspolrzednych [1-wysokosc] -> [0-1]
jj = (j-1)/(iy-1);
% B^y_j(y)
funy_tab(bj,j) = compute_spline(knot_vectory,py,bj,jj);
endfor
endfor
preprocess=toc
tic;
% petla po funkcjach
for bi = 1:nx
for bj = 1:ny
% petla po pixelach obrazka %na których funkcje sa niezerowe
for i=xx(knot_vectorx(bi)):xx(knot_vectorx(bi+px+1))
for j=yy(knot_vectory(bj)):yy(knot_vectory(bj+py+1))
% B^x_i(x)
funi = funx_tab(bi,i);%compute_spline(knot_vectorx,px,bi,ii);
% B^y_j(y)
funj = funy_tab(bj,j);%compute_spline(knot_vectory,py,bj,jj);
ff = funi*funj;
R1(i,j) = R1(i,j) + floor(ff*RR((bj-1)*nx+bi));
G1(i,j) = G1(i,j) + floor(ff*GG((bj-1)*nx+bi));
B1(i,j) = B1(i,j) + floor(ff*BB((bj-1)*nx+bi));
endfor
endfor
endfor
endfor
RGB=X;
RGB(:,:,1) = R1;
RGB(:,:,2) = G1;
RGB(:,:,3) = B1;
rebuild=toc
imshow(RGB);
function resx=xx(x)
resx = floor((ix-1)*x+1);
endfunction
function resy=yy(y)
resy = floor((iy-1)*y+1);
endfunction
function val=bitmp(M,x,y)
val = zeros(size(x));
for i=1:size(x,1)
for j=1:size(x,1)
val(i,j)=M(xx(x(1,i)),yy(y(1,j)));
endfor
endfor
endfunction
%funkcja wyliczajaca stopien wielomianow
function p=compute_p(knot_vector)
%pierwszy wpis w knot_vector
initial = knot_vector(1);
%dlugosc knot_vector
kvsize = size(knot_vector,2);
p = 0;
%sprawdzenie ilosci powtorzen pierwszego wpisu w knot_vector
while (p+2 <= kvsize) && (initial == knot_vector(p+2))
p = p+1;
endwhile
return
endfunction
%funkcja sprawdzajaca poprawnosc knot_vector
function t=check_sanity(knot_vector,p)
initial = knot_vector(1);
kvsize = size(knot_vector,2);
t = true;
counter = 1;
%jesli ilosc powtorzen na poczatku nie jest zgodna ze stopniem wielomianow
for i=1:p+1
if (initial != knot_vector(i))
%zwroc falsz
t = false;
return
endif
endfor
%jesli w srodku knot_vector jest za duzo powtorzen
for i=p+2:kvsize-p-1
if (initial == knot_vector(i))
counter = counter + 1;
if (counter > p)
%zwroc falsz
t = false;
return
endif
else
initial = knot_vector(i);
counter = 1;
endif
endfor
initial = knot_vector(kvsize);
%jesli ilosc powtorzen na poczatku nie jest zgodna ze stopniem wielomianow
for i=kvsize-p:kvsize
if (initial != knot_vector(i))
%zwroc falsz
t = false;
return
endif
endfor
%jesli kolejny element knot_vector jest mniejszy niz poprzedni
for i=1:kvsize-1
if (knot_vector(i)>knot_vector(i+1))
%zwroc falsz
t = false;
return
endif
endfor
return
endfunction
%funkcja wyliczajaca rekurencyjnie funkcje bazowa zgodnie z wzorem Cox-de-Boor
function y=compute_spline(knot_vector,p,nr,x)
%funkcja (x-x_i)/(x_{i-p}-x_i)
fC= @(x,a,b) (x-a)/(b-a);
%funkcja (x_{i+p+1}-x)/(x_{i+p+1}-x_{i+1})
fD= @(x,c,d) (d-x)/(d-c);
%x_i
a = knot_vector(nr);
%x_{i-p}
b = knot_vector(nr+p);
%x_{i+1}
c = knot_vector(nr+1);
%x_{i+p+1}
d = knot_vector(nr+p+1);
%funkcja liniowa dla p=0
if (p==0)
y = 0 .* (x < a) + 1 .* (a <= x & x <= d) .+ 0 .* (x > d);
return
endif
%B_{i,p-1}
lp = compute_spline(knot_vector,p-1,nr,x);
%B_{i+1,p-1}
rp = compute_spline(knot_vector,p-1,nr+1,x);
%(x-x_i)/(x_{i-p)-x_i)*B_{i,p-1}
if (a==b)
%jesli wezly w knot_vector sie powtarzaja trzeba to uwzglednic
y1 = 0 .* (x < a) + 1 .* (a <= x & x <= b) + 0 .* (x > b);
else
y1 = 0 .* (x < a) + fC(x,a,b) .* (a <= x & x <= b) + 0 .* (x > b);
endif
%(x_{i+p+1}-x)/(x_{i+p+1)-x_{i+1})*B_{i+1,p-1}
if (c==d)
%jesli wezly w knot_vector sie powtarzaja trzeba to uwzglednic
y2 = 0 .* (x < c) + 1 .* (c < x & x <= d) + 0 .* (d < x);
else
y2 = 0 .* (x < c) + fD(x,c,d) .* (c < x & x <= d) + 0 .* (d < x);
endif
y = lp .* y1 .+ rp .* y2;
return
endfunction
% Ilosc elementow w knot vector
function n = number_of_elements(knot_vector,p)
initial = knot_vector(1);
kvsize = size(knot_vector,2);
n = 0;
for i=1:kvsize-1
if (knot_vector(i) != initial)
initial = knot_vector(i);
n = n+1;
endif
endfor
endfunction
% Tworzy prosty knot vector bez powotrzen w srodku
function knot = simple_knot(elems, p)
pad = ones(1, p);
knot = [0 * pad, 0:elems, elems * pad];
knot = knot/elems;
endfunction
% Ilosc funkcji bazowych w knot wektorze
function n = number_of_dofs(knot,p)
n = length(knot) - p - 1;
endfunction
function first = first_dof_on_element(knot_vector,p,elem_number)
[l,h] = element_boundary(knot_vector,p,elem_number);
first = lookup(knot_vector, l) - p;
endfunction
function [low,high] = element_boundary(knot_vector,p,elem_number)
initial = knot_vector(1);
kvsize = size(knot_vector,2);
k = 0;
low=0;
high=0;
for i=1:kvsize
if (knot_vector(i) != initial)
initial = knot_vector(i);
k = k+1;
endif
if (k == elem_number)
low = knot_vector(i-1);
high = knot_vector(i);
return;
endif
endfor
endfunction
% Zwraca zakres (indeksy) funckji bedacych niezerowymi na zadanym wektorze wezlow
function [low,high] = dofs_on_element(knot_vector,p,elem_number)
low = first_dof_on_element(knot_vector,p,elem_number);
%poniewaz mamy dokladnie p+1 niezerowych funkcji nad elementem
high = low + p;
endfunction
% Row vector of points of the k-point Gaussian quadrature on [a, b]
function xs = quad_points(a, b, k)
% mapowanie punktow
map = @(x) 0.5 * (a * (1 - x) + b * (x + 1));
switch (k)
case 1
xs = [0];
case 2
xs = [-1/sqrt(3), ...
1/sqrt(3)];
case 3
xs = [-sqrt(3/5), ...
0, ...
sqrt(3/5)];
case 4
xs = [-sqrt((3+2*sqrt(6/5))/7), ...
sqrt((3-2*sqrt(6/5))/7), ...
sqrt((3-2*sqrt(6/5))/7), ...
sqrt((3+2*sqrt(6/5))/7)];
case 5
xs = [-1/3*sqrt(5+2*sqrt(10/7)), ...
-1/3*sqrt(5-2*sqrt(10/7)), ...
0, ...
1/3*sqrt(5-2*sqrt(10/7)), ...
1/3*sqrt(5+2*sqrt(10/7))];
otherwise
xs = [-1/3*sqrt(5+2*sqrt(10/7)), ...
-1/3*sqrt(5-2*sqrt(10/7)), ...
0, ...
1/3*sqrt(5-2*sqrt(10/7)), ...
1/3*sqrt(5+2*sqrt(10/7))];
endswitch
xs = map(xs);
endfunction
% Row vector of weights of the k-point Gaussian quadrature on [a, b]
function ws = quad_weights(a, b, k)
switch (k)
case 1
ws = [2];
case 2
ws = [1, 1];
case 3
ws = [5/9, ...
8/9, ...
5/9];
case 4
ws = [(18-sqrt(30))/36, ...
(18+sqrt(30))/36, ...
(18+sqrt(30))/36, ...
(18-sqrt(30))/36];
case 5
ws = [(322-13.0*sqrt(70))/900, ...
(322+13.0*sqrt(70))/900, ...
128/225, ...
(322+13.0*sqrt(70))/900, ...
(322-13.0*sqrt(70))/900];
otherwise
ws = [(322-13.0*sqrt(70))/900, ...
(322+13.0*sqrt(70))/900, ...
128/225, ...
(322+13.0*sqrt(70))/900, ...
(322-13.0*sqrt(70))/900];
endswitch
% Gaussian quadrature is defined on [-1, 1], we use [a, b]
%ws = ws * (b-a)/2;
endfunction
endfunction
Autorzy kodów w MATLABie: Marcin Łoś i Maciej Woźniak.